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Abstract Homogeneous and isotropic models are studied in the Jordan 
frame of the second order gravity theory. The late time evolution of the 
models is analysed with the methods of the dynamical systems. The normal 
form of the dynamical system has periodic solutions for a large set of initial 
conditions. This implies that an initially expanding closed isotropic universe 
may exhibit oscillatory behaviour. 

Keywords Isotropic cosmologies Higher-order gravity Dynamical systems 



1 Introduction 

Quadratic gravitational Lagrangians were proposed shortly after the formu- 
lation of general relativity (GR) as alternatives to Einstein's theory. Gravity 
modifications in the form of higher-order curvature invariants in the La- 
grangian are generally known as higher-order gravity (HOG) theories. They 
arise in string-theoretic considerations, e.g., brane models with Gauss-Bonett 
terms Q] or models with a scalar field coupled to the Gauss-Bonett invari- 
ant [2] (see [3] for a review) and generally involve linear combinations of all 
possible second order invariants that can be formed from the Riemann, Ricci 
and scalar curvatures. A quarter of a century ago there was a resurgence 
of interest in such theories in an effort to explain inflation. The reasons for 
considering HOG theories were multiple. Firstly, it was hoped that higher 
order Lagrangians would create a first approximation to quantum gravity, 
due to their better renormalisation properties than GR [4]. Secondly, it was 
reasonable to expect that on approach to a spacetime singularity, curvature 
invariants of all orders ought to play an important dynamical role. Far from 
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the singularity, when higher order corrections become negligible, one should 
recover GR. Furthermore, it was hoped that these generalized theories of 
gravity might exhibit better behavior near singularities. Thirdly, inflation 
emerges in these theories in a most direct way. In one of the first inflation- 
ary models, proposed in 1980 by Starobinsky [5], inflation is due to the R 2 
correction term in a gravitational Lagrangian L = R + pR 2 where j3 is a 
constant. 

Recently there is a revival of interest in HOG theories in an effort to ex- 
plain the accelerating expansion of the Universe J6j[7] . The general idea is to 
add an 1/ R term to the Einstein-Hubert Lagrangian or more generally to con- 
sider R+aR~ n Lagrangians [HUH]- As the Universe expands, one expects that 
the inverse curvature terms will dominate and produce the late time acceler- 
ating expansion. Most studies are restricted to simple Friedmann-Robertson- 
Walker (FRW) models because of the complexity of the field equations. At 
present, the observational viability of these models is a subject of active re- 
search (see [IO l [TT | ri2 j ri3 ] and references therein). However, it seems that a 
large class of Lagrangians may fit the observational data, but simple models 
based on the FRW metric are insufficient to pick the correct Lagrangian (see 
for example [T3] for the reconstruction of the f(R) theory which best repro- 
duces the observed cosmological data). For more general spacctimes it may 
even be meaningless to say that R is small in some epoch of cosmic evolution 
and large in some other one (for a thorough critic see [15]). 

In this paper we investigate the late time evolution of flat and positively 
curved FRW models with a perfect fluid in the R + (3R 2 theory. This is the 
simplest generalization of the Einstein-Hilbert Lagrangian and the addition 
of the quadratic term represents a correction to general relativity. The simple 
vacuum case was studied in [IB] , where oscillatory behaviour of the solutions 
of closed models was found. Since HOG theories in vacuum are conformally 
equivalent to GR with a scalar field, it is tempting to say that the R 2 con- 
tribution has predictable cosmological consequences [IT] . However, this is an 
oversimplification of the picture (see [TUlfT5] for specific examples). The two 
frames are mathematically equivalent, but physically they provide different 
theories. In the Jordan frame, gravity is described entirely by the metric 
g^. In the Einstein frame, the scalar field exhibits a non- metric aspect of 
the gravitational interaction, reflecting the additional degree of freedom due 
to the higher order of the field equations in the Jordan frame. Inclusion of 
additional matter fields, further complicates the situation and while the field 
equations in the Einstein frame are formally the Einstein equations, never- 
theless this theory is not physically equivalent to GR. There is no universally 
acceptable answer to the issue "which conformal frame is physical" [19] (see 
[20] for a thorough analysis of different views). 

The plan of the paper is as follows. Next Section contains a short com- 
ment on the stability of well-known power-law solutions. The field equations 
are written as a constrained four-dimensional polynomial dynamical system. 
Section 3 contains the analysis of the flat case. The so-called normal form 
of the dynamical system greatly simplifies the problem, since two of the 
equations decouple. In Section 4 we study the qualitative behaviour of the 
solutions near the equilibrium points of positively curved models and analyse 



3 



their late time evolution. It is shown that an initially expanding closed FRW 
universe may exhibit oscillatory behaviour. 

2 Field equations 

The general gravitational Lagrangian in four-dimensional spacetimes con- 
tains curvature invariants of all orders, R, R 2 , R^R^ , .... The term 
Rfj.vpaR^'" 7 is omitted because of the Gauss-Bonnet theorem. A further sim- 
plification can be done in homogeneous and isotropic spacetimes where the 
variation of R^R^ with respect to the metric is proportional to the varia- 
tion of R 2 [21]. We conclude that for isotropic cosmologies the gravitational 
Lagrangian contains only powers of the scalar curvature and we may consider 
HOG theories derived from Lagrangians of the form 

L = f (R) V-9 + ^matter, 

where / is an arbitrary smooth function. It is well-known that the corre- 
sponding field equations are fourth-order and take the form 

/' (R) R pv - i/ (R) 9ltu - V M V„/' (R) + g^Df (R) = T^, (1) 

where □ = g a PV a Vp and a prime (') denotes differentiation with respect to 
R. The generalised Bianchi identities imply that V^T^ V = 0. Contraction of 
(P) yields the trace equation 

3Df (R) + f (R) R -2f (R) = T. (2) 

In contrast to GR where the relation of R to T is algebraic, in HOG theories 
the trace equation ([2]) is a differential equation for R, with T as source term 
[22] . This suggests that in HOG theories both the metric and the scalar 
curvature are dynamical fields. 

In the following, we consider a quadratic Lagrangian without cosmolog- 
ical constant, i.e. f (R) = R + PR 2 , (3 > 0, and confine our attention to 
cosmologies with a perfect fluid with energy density p and pressure p, of the 
form 

p=(7-l)p, 0< 7 <2. 

For homogeneous and isotropic spacetime^] described by the standard FRW 
metric we need the following useful relations = 1, 2, 3 and k = 0, ±1) 

a „ (a a 2 k\ „ /da 2 k \ 

R 0Q = -3-, R l3 = - + 2— + 2— ) gij , R : = 6 - + — + — . 3 

The 00 component of (Q]) is 

(4) 



1 We adopt the metric and curvature conventions of [55]. Here, a (t) is the scale 
factor, an overdot denotes differentiation with respect to time t, and units have 
been chosen so that c = 1 = 8nG. 
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where H = a/a, is the Hubble function. 

At this point we make a digression. For flat, k = 0, models, differentiating 
the relation R = 6H +12H 2 (which comes from the third of ([3])) with respect 
to t, equation (HJ) takes the form 

H 2 + 6/3 (2HH + 6H 2 H - H 2 ^j = ^p. (5) 

For radiation, 7 = 4/3, there exists a special solution a (t) = t 1 / 2 as t — > 0. As 
Barrow and Middleton point out [24] , this is also an exact vacuum solution of 
the purely quadratic theory, in the sense that it solves 2HH+6H 2 H—H 2 = 0. 
However, this solution is unstable as we shall see in a moment (see |25] for 
a detailed stability analysis of isotropic models in general / (R) theories). 
Following Carroll et al [8], we reduce the order of ([5]) in vacuum by defining 

X = -H, Y = H => H = -Y^-. 

aJC 

Then, the asymptotic values of the function U (X) = —X 2 /Y as X — > 0, 
correspond to the exponents p for power-law solutions a (t) = t p . We apply 
this technique to §E§ and we find the first-order equation 

dX = YTpX^- i x\ lJ -2)^ (6) 

with the corresponding direction field shown in Fig. 1. We see that the so- 
lution U (X) = 1/2 is a past attractor (X — * —00), but becomes unstable 
as X — > 0, in agreement with [25], and that |{7(X)| — ► 0, corresponding to 

the singularity H — * 00. This is also evident by studying the asymptotic 

behaviour of solutions of the linearised equation near the constant solutions 
1/2 and 0. 

We now continue our discussion about the choice of variables. Setting 
x = X/a, the 00 equation becomes 



H z + kx z + 2(3 



R 2 

R (H 2 + kx 2 ) + HR - — 



and the evolution equation for x is 

x = -xH. (8) 
The evolution equation for H comes from the third of ^ and takes the form 

H = -R- 2H 2 - kx 2 . (9) 
6 

The conservation equation 

p = -3 1P H, (10) 

is a consequence of the Bianchi identities. With the relation OR = —R—3HR, 
equation © becomes 

i?. + 3tf7?+^i?=^(4-3 7 )p. (11) 
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Fig. 1 Direction field of the reduced equation (0 



This equation is usually considered as superfluous, since it follows from the 
differentiation of the equation ([7]) with respect to t; for example equation 
(fTTjl was used in [26] as a control of the accuracy of numerical investigations 
of Bianchi type I and IX models. Thus, one can chose R,x,H and p as 
dynamical variables which obey the evolution equations 0, ®, © and 
(fTtJf and constitute a four-dimensional dynamical system. 

The form of equation ([TJ suggest the choice of expansion normalized vari- 
ables of the type 

u~R/H, f2~ p/H 2 ,... (12) 

Detailed studies using this approach for general / (R) theories can be found 
in 12,27 . However, even in the simplest case of flat models in vacuum, 
numerical investigation of the system with initial values H (0) > shows 
that H (t) exhibit damped oscillations, with almost zero minima. Therefore 
although permissible, the transformation (|12p may induce fake singularities 
to the solutions. Moreover, as mentioned in [25] . a drawback of this choice of 
variables is that it does not give a complete description of the evolution for 
bouncing or recollapsing models. If the Hubble parameter passes through zero 
the logarithmic time coordinate is ill-defined and the transformation (fT2|) is 
singular (see [28^ for a comparison of compact and non-compact variables). 

In order to circumvent these difficulties, we introduce one more degree 
of freedom by using (jlip . thus augmenting the dimension of the dynamical 
system. The state (R, R, x, p, H) of the system lies on the hypersurface of 
M 5 defined by the constraint ([7]) . The presence of (R, R) in the state vector 
reflects the fact that there are additional degrees of freedom in HOG theories 
than in GR (cf. the remark after equation ([2])). 
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We define x\ — R, x 2 = R and use ((7]) to eliminate p, so that equations 
(fTTjl , JHJ and ([9|) constitute a four-dimensional system. The parameter f3 may 
be used to define dimcnsionless variables by rescaling 



S!-^!, ^^V^ X2 ' ^7P' ^V^' *- > ^ t ' 

and our system becomes 

±1 = x 2 , 

z 2 = -H + (1 - 3 7 ) x 2 H + (4 ~ 37) Z, (13) 
i = —xH, 

H = 2xi - 2H 2 - kx 2 , 

with Z = [H 2 + kx 2 - + 4xi (H 2 + kx 2 )] . 

Remark. The system (fT3"|) is not an arbitrary "free" four-dimensional 
system. In view of ([7]) the initial conditions have to satisfy the condition 

H 2 + kx 2 + Axi (H 2 + kx 2 ) + AHx 2 - Ax\ > 0. (14) 



With a little manipulation of the equations (|T3[) it can be shown that, once 
we start with initial conditions satisfying (fl"4"]) at time to , the solutions of the 
system satisfy this inequality for all t > t - Thus the field equations share 
the general property of the Einstein equations, namely that the subsequent 
evolution of the system is such that the solutions respect the constraint. 



3 Flat models 

In the flat, (k — 0), case the dimension of the dynamical system is reduced 
by one, since the evolution equation for x decouples from the remaining 
equations. The corresponding system is 



2 , (4-37) „ 2 . , u2 



xi = x 2 , 

x 2 = -xi + (1 - 3 7 ) x 2 H - (4 - 37) x'i + K - ^" H 2 + (4 - 3 7 ) Xl H 2 

H = 2xi - 2H 2 , (15) 

i.e., the vector field does not depend on the dynamical variable x. Vacuum 
models are significantly simpler to analyse, but we do not study them sep- 
arately as they arise formally by setting 7 = 4/3 in all the equations, while 
(I14p describing the phase space becomes equality. The only equilibrium point 
of (|15p is the origin and corresponds to flat empty spacetime. The eigenval- 
ues of the Jacobian matrix at the origin are ±i,0 and therefore we cannot 
infer about stability using the linearisation theorem. For nonhyperbolic equi- 
librium points there exist no general methods for studying their stability. 
The normal form of the system may provide some information about the 
behaviour of the solutions near the equilibrium. The normal form theory 
consists in a nonlinear coordinate transformation that allows to simplify the 
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nonlinear part of the system (cf. [29] for a brief introduction). This task will 
be accomplished in three steps in some detail for the convenience of readers 
with no previous knowledge of the method. 

1. Let P be the matrix formed from the eigenvectors which transforms 
the linear part of the vector field into Jordan canonical form. We write (|15p 
in vector notation (with x = (x\, x 2 , H)) as 



x = Ax + F (x) , 



(16) 



where A is the linear part of the vector field and F (0) = 0. 

2. Using the matrix P, we define new variables, (yi,y 2l y) = y, by the 
equations 

Vi=xi, y 2 = -x 2 , y = H + 2x 2l 
or in vector notation y = fx, so that (fTB")) becomes 



y = p- x APy + P _1 F (Py) . 
Denoting the canonical form of A by B wc finally obtain the system 

y = Py + f(y), (17) 
where f (y) := P _1 F (Py) . In components system (fl~7f is 



m 




0-10 




yi 


m 




1 




2/2 


y 









y 





(4 - 3 7 ) yl - (37 + 2) y\ - ^y 2 - 3yy 2 - (4 - 3 7 ) yi (2y 2 + yf 
-2 (4 - 37) yi + 2 (3 7 - 2) y\ - f y 2 - 2yy 2 + 2 (4 - 3 7 ) Vl (2y 2 + yf 



Inequality (|14p imposes the constraint 

y 2 - 4 (y 2 + y 2 ) + A Vl {y + 2y 2 f > 0. 
3. Under the non-linear change of variables 



(18) 



m + 3 7 y? + (37 - 2) y 2 + 



4-3 7 2 

— a — y 



4 y2V > 



2/2 + 4ym + jyiy, 



1)2 

V^y- 22/iJ/2 + 22/i2/, 
and keeping only terms up to second order, the system transforms to 



(19) 



m = -2/2- -2/ij/ + O (3) , 



2/2=2/1- 2^22/ + 0(3), 
2/ = 6(7-l)(2/?+2/2 2 )- |W<3(3). 
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Finally, defining cylindrical coordinates (y\ = r cos 9, y 2 = rsui9,y — y) , we 
obtain 

f = ~ry + 0(3), 

9 = 1 + O (2) , (20) 
y = 6( 7 -l)r 2 -|V + 0(3). 
We may continue to simplify the third order terms and the result should be 

f = airy + a 2 r 3 + a 3 ry 2 + O (4) , 

= 1 + 0(2), 

1 = hr 2 + b 2 y 2 + b 3 r 2 y + b 4 y 3 + 0(4). 

This is the normal form in cylindrical coordinates of every three-dimensional 
vector field with linear part 



"0 -1 o" 




2/i 


1 




2/2 







V 



(see [3D] P- 377). However, since we are interested on the behaviour of the 
solutions only near the origin, we truncate the vector field at O (2) . We note 
that the 9 dependence of the vector field has been eliminated, so that we 
can study the system in the (r, y) space. The second equation of (|2H)l implies 
that the trajectory in the y\ — y 2 plane spirals with angular velocity 1. The 
projection of (|20|) on the r — y plane is 

3 

r = ~2 ry > 

y = 6( 7 -l)r 2 -|V. (21) 

This system belongs to a family of systems studied in 1974 by Takens [31] 
(see also [30 for a description of all phase portraits for different values of the 
parameters). 

System (|2ip is invariant under the transformation t — > —t, y — * — y (which 
implies that all trajectories are symmetric with respect to the r axis) and 
the line r = is invariant. Note also that the system ([2T|) has invariant 
lines y — ±2r (compare with the similar system in |32j). The behaviour of 
the solutions depends on the parameter 7 and as we shall see, 7 = 1 is a 
bifurcation value. 

Case I, 7 < 1. We observe that y is always decreasing along the orbits 
while r is decreasing in the first quadrant. Since no trajectory can cross 
the line y = 2r, all trajectories starting above this line, approach the origin 
asymptotically. The phase space of the dynamical system (|2"Tj) is not the 
whole r — y plane, because of the constraint (|18|) . In terms of the variables 
(fl~9)) and neglecting fourth-order terms the constraint becomes 



(y 2 -4r 2 )(l + 6 7 yi)>0. 



(22) 
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Fig. 2 Phase portrait of (|21[) for 7 < 1 and 7 > 1. The invariant lines y — ±2r 
are shown with heavy lines and y = 2V7 — 1 r by the dotted line. 



Since we are interested on trajectories starting close to the origin, for initial 
values of y\ satisfying |j/i (0)| < I/67 the first of ([21]) guarantees that \y\ (t)\ 
remains less than I/67, so that (|22[) implies 

y 2 - 4r 2 > 0. 

Therefore, we should consider only trajectories starting above the line y = 2r 
and according to the previous discussion all these trajectories asymptotically 
approach the origin. 

Case II, 7 > 1. In the first quadrant r is decreasing along the orbits and, 
that y vanishes along the line y = 2^/7 — 1 r. Once a trajectory crosses the 
line y = — 1 r, it is trapped between the lines y — 2y/j — 1 r and y = 2r, 
and since r < 0, it approaches the origin asymptotically. 

Case III, 7 = 1. It is evident that all trajectories are straight lines ap- 
proaching asymptotically the origin. 

We conclude that the late time behaviour of flat models is similar to the 
future predicted by GR. More precisely, all initially expanding flat models 
close to the state R = R = H = p = 0, asymptotically approach the flat 
empty spacetime. 



4 Positively curved models 

In this section we consider an initially expanding closed universe described 
by the full four-dimensional system (|I3p with k = I. There are two equilibria: 
M. : xi = 0, x 2 — 0, x = 0, H — 0. This corresponds to the limiting 
state of an almost empty, slowly varying universe with H — * 0, while the 
scale factor goes to infinity. The point M. which resembles to the Minkowski 
solution, is located at the boundary of the phase space. 
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S : H = 0, x 2 = 0, x = x = yj^=i,xi = #i = x 2 /2, with 2/3 < 7 < 4/3. 
It is a static solution and the eigenvalues of the Jacobian matrix at S are 



8 - 9 7 ± y/3-y (36 7 2 - 69 7 + 32) 
2 (3 7 ~ 4) ■ 

The real parts of the eigenvalues are nonzero for almost alQ permissible 
values of 7, and we conclude that the local stable and unstable manifolds 
through S arc both two-dimensional. The point S corresponds to the Einstein 
static universe, where the effective cosmological constant is provided by the 
curvature equilibrium x\. To see this, it is sufficient to write equation (0 in 
the original variables at the equilibrium point as 

kx + - = -p, 

with A > 0. The cosmological constant depends on both the parameters (3 
and 7, for example, for 7 = 1, A = (12/3— l) 2 /2f3. Static solutions have 
little interest as cosmological models and we turn our attention to the other 
equilibrium, M.. 

The point M. is a nonhyperbolic equilibrium and we find again the normal 
form of the system, which is given by (|A.29|) in the Appendix. Defining 
cylindrical coordinates — r cos 0, y 2 = r sin 9,x = x,y = y) , we obtain 

r = --ry + 0(3), 
6 = 1 + 0(2), 

x = -xy + O (3) , (23) 
^6( 7 -l)r 2 -^, 2 -|V + 0(3). 

We truncate the vector field at O (2) and we note again that the 8 dependence 
of the vector field has been eliminated, so that we can study the system in the 
(r, x, y) space. We write the first and third of (|23|) as a differential equation 

dr 3 r 
dx 2 x ' 

which has the general solution 

r = Ax 3/2 , A>0. (24) 

2 More precisely, the real parts of the eigenvalues are zero for 

2 ^ 23 - vTf . „ n 
3<7<^^-0.79, 

and nonzero in the rest of the interval (|, |). 
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Substitution of (f24|) into the fourth equation of ([23]) yields the projection of 
the fourth-dimensional system on the x — y plane, namely 

x = -xy, 

y = 6( 7 -l)x 3 -^x 2 -^y 2 , b > 0. (25) 

Some general properties of the solutions of ([^5)1 follow by inspection. By 
standard arguments all trajectories are symmetric with respect to the x axis. 
Note that the new x defined by (|A.28[) remains non-negative for initial values 
2/i (0) sufficiently small and, since the line x = is invariant, any trajectory 
starting at the half plane x > remains there for all t > 0. System (|25j) has 
two equilibrium points, the origin (0,0) and (a;*,0), where 

3 7 - 2 
X *~ 26(7-1)' 

and therefore x* > for 7 < 2/3 or 7 > 1. The origin is again a nonhyperbolic 
equilibrium point. Computation of the Jacobian matrix at the equilibrium 
(x*,0) shows that for the linearised system, this point is a saddle for < 
7 < 2/3, and a center for 1 < 7 < 2. For 7 > 1, it is easy to see that x is 
decreasing in the first quadrant and y is decreasing along the orbits in the 
strip < x < x*. On any orbit starting in the first quadrant with x < x», y 
becomes zero at some time and the trajectory crosses vertically the a;— axis. 
Once the trajectory enters the second quadrant, x increases and y decreases. 
For 2/3 < 7 < 1, all trajectories starting in the first quadrant follow the 
same pattern. 

System (|25[) has a first integral, viz. 

< P {Xl y) = - 2 ± X ^) +X ^ + ^- 1 7^1, 

<P(x,y) = - + ^3, 7 = !• 

X x° 

This can be seen by writing (|25]) as 

dy _ bx 3 - ^x 2 - %y 2 
dx —yx 

Setting y 2 — z, we obtain a linear differential equation for z which is easily 
integrable. The level curves of 4> are the trajectories of the system. 

We shall show for the system ([^5)1 that: (i) {or every 7 e (|, 2] there are 
no solutions asymptotically approaching the origin (ii) for 7 G (§>!] there 
are no periodic solutions and (Hi) for 76 (1,2] there exist periodic solutions 
and the basin of every periodic trajectory is the set 

l (x, y) £ R 2 : y 2 + x 2 - y x 3 < 0, x > j . 
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Proof The proof mimics that found in [16j for the simple case 7 = 4/3. 
Let (p (x, y) = C We have 

f = * l (-l + jX + C a *>-*y (26) 

which implies that the function 

must be non-negative. We consider two cases. 

1. C > 0. Then / is strictly increasing for x > and / (0) = — 1, thus / 
has a unique root 27 > depending on C. It follows that for C > any orbit 
starting in the first quadrant satisfies 

X > X! (C) > 0, 

i.e., there are no solutions approaching the axis x — 0. These solutions are 
not closed since they intersect the x— axis only once at 27 (C). 

2. C < 0. If 2/3 < 7 < 1, then / has again a unique root hence the trajec- 
tories follow the same pattern as in case 1. If 7 > 1, then / has two zeros, say 
X\{C) < x%(C). This means that < 27(C) < x < x 2 {C), i.e., x is bounded, 
and by (f2"6"|) , so is y. Thus, an orbit of (|2"5j) starting in the first quadrant 
crosses the x— axis at x± (C) and re-enters in the first quadrant crossing the 
x— axis at xi (C) i.e. it is a closed curve and represents a periodic solution. 
The curve corresponding to C — separates the phase space into two dis- 
joint regions I and II. In region II, (C < 0), every trajectory corresponds to a 
periodic solution and we conclude that the basin of every periodic trajectory 
is the set y 2 + x 2 - |a; 3 < 0. 

Remark. The mere existence of closed orbits around the equilibrium 
point (x*, 0) could be inferred from the following theorem: If an equilibrium 
point p is a center for the linearised system and all trajectories are symmetric 
with respect to the x axis, then p is also a center for the nonlinear system 
(f25|) (cf. [29] Theorem 6, page 141). In the above proof we also determine the 
subset of the phase space which contains all periodic orbits. 

Using all this information we may sketch the phase portrait of the system 
(Fig. 3). 

For 7 < 2/3 the homoclinic curves to the origin as well as the curves 
approaching the origin indicate that an initially expanding closed universe 
may avoid recollapse. This result is also valid in GR when matter fields violate 
the strong energy condition [33 . For 7 £ (|, 1] every solution curve becomes 
unbounded and we may interpret this as an indication that this universe 
recollapses. 

The range (1,2] for 7 is the more interesting because of the periodic orbits 
in region II. The phase portrait in Fig. 3 may lead to the conclusion that the 
periodic orbits are far from the origin. However, the position of the cycles 
in the phase space depends on the constants b in ([2"5]) and C in (|2"6"j) and 
therefore, for suitable values of b and C, there exist periodic orbits arbitrary 
close to the origin. Note that the periodic solutions of ([25| induce periodicity 
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to the full four-dimensional system (|23[) or (|A.29|) . In fact, if x (t) and y (t) 
are periodic solutions, then (f2"4"| implies that the solutions 

2/1 (t) = r (t) cos (t + 0„ ) , 2/2 (t) = r (t) sin (t + 0„) 

oscillate in the j/i — y% plane with a periodic amplitude r(t). Note also that 
the periodic motion in the x — y plane is independent from the rotation of r 
in the j/i — J/2 plane. 

Obviously we cannot assign a physical meaning to the new variables 
(yi,y2,x,y) , since the transformations (|A.27[) and (|A.28[) have "mixed" the 
original variables of (|13[) in a nontrivial way. However, the periodic character 
of the solutions of (IA.29j) whatever the physical meaning of the variables be, 
has the following interpretation. Close to the equilibrium of the original sys- 
tem H3\) , there exist periodic solutions. This result is in agreement with [34] 
where it is shown that for Lagrangians R + (3R rn , bounces of closed models 
are allowed for every integer value of m. 

Remark. As mentioned in Section 2 numerical experiments show that 
the solutions of the system have oscillatory behavior. This property is intu- 
itively evident by looking at the harmonic oscillator, equation Using the 
rescalings of the dynamical variables along with the scaling p —f [3p 1 equation 
(|TT|) becomes 

R + 3HR + R= ( 4 ~ 3 ^) p; 

which is the equation of a forced, damped harmonic oscillator with unit 
angular frequency. Qualitative arguments supported by numerical solutions 
indicate that oscillatory motion, possibly slightly damped, is essentially in- 
dependent of k and 7 and becomes the late time behaviour. This is revealed 
in the normal forms of the systems, ([2H)l and (f2"3")) . where the equation 
shows oscillatory motion with unit angular frequency. We emphasize again 
that this motion is independent from the periodic motion in the x—y plane. 
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The normal form analysis reveals both kinds of oscillatory behaviour and dis- 
tinguishes the models which actually exhibit undamped oscillations at late 
time. 

5 Discussion 

We analysed the qualitative behaviour of flat and positively curved FRW 
models filled with ordinary matter described by a perfect fluid in the Jordan 
frame of the R + (3R 2 theory. We have shown that initially expanding flat 
models close to the equilibrium H — 0, i? = 0, p = are ever expanding 
and asymptotically approach flat empty spacetime. Therefore, the late time 
behaviour of flat models is a common property of quadratic gravity and 
GR. Closed models can avoid recollapse for 7 < 2/3, but not for 7 in the 
range [|, 1] , thereby also behaving similarly to the corresponding general- 
relativistic cosmologies. 

The interesting feature is the existence of periodic solutions near the ori- 
gin for 7 > 1. This is not revealed in the Einstein frame (see [32]), possibly 
because in that investigation, the scalar field related to the conformal trans- 
formation is not coupled to matter, i.e., the matter Lagrangian is added 
after performing the conformal transformation. Were the two frames phys- 
ically equivalent then, a naive physical explanation of the cyclic behavior 
could rely on the scalar field which behaves like a "cosmological constant" in 
the high curvature limit. The perfect fluid which is also present dominates 
in the low curvature regime allowing for a recollapse, but then the effective 
cosmological constant induces a bounce in the high-curvature regime. 

The periodic solutions imply that an initially expanding closed universe 
can avoid recollapse through an infinite sequence of successive expansions 
and contractions. The oscillatory open model proposed by Stcinhardt and 
Turok |35j has renewed interest in cyclic universes. However, observations 
do not exclude ^totai to be slightly larger than one [3B]. Oscillatory closed 
models were considered in the context of Loop Quantum Cosmology [37] and 
oscillatory (but not periodic) solutions in GR were found in 38] for closed 
models containing radiation and dust or scalar field. In Fig. 3 the basin of all 
periodic trajectories of the ([23]) is the domain on the right of the C — curve 
and since it is an open subset of the phase space, we conclude that there is 
enough room in the set of initial data of (fl"3]) which lead to an oscillating scale 
factor. The R + f3R 2 theory has offered a successful inflationary model, but 
whether it is capable to explain the acceleration of the universe is an open 
question that needs to be studied in more detail. In particular, the observed 
slow acceleration must be related to the periods of the closed curves of ([25]) . 

The normal form theory is a powerful method for determining the qual- 
itative behaviour of a dynamical system near a nonhyperbolic equilibrium 
point, but does not give any information about the structure of the solu- 
tions far from this equilibrium. Our results are based on an analysis of the 
behaviour of the dynamical system (|13[) only near the equilibrium solutions. 
The geometry of the trajectories of a four-dimensional dynamical system may 
be quite complicated, e.g. strange attractors may be present. For the system 
(|13p the whole picture may come in view only from the investigation of the 
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global structure of the solutions. The study of this question is an interesting 
challenge for mathematical relativity. 
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Appendix: Normal form of (113[) 



Following the usual algorithm, we use the matrix which transforms the linear 
part of (|13p into Jordan canonical form and define new variables by 

yi=xt, V2 = -X2, x = x, y = H + 2x 2 , (A. 27) 

so that (fl"3l) becomes 



~yx~ 




"0 -1 0" 




~yi~ 




m 




1 00 




2/2 




X 









X 


+ 


. V . 









. y - 






(4- 


37)y?-(37 


+ 2)y 2 2 - 


-2(4- 


3 7 )2/ 2 + 2 (3 7 - 


2)2/1 



c 2 + y 2 ) - '3yy 2 - (4 - 3 7 ) Z 3 



-2xy 2 - xy 



with Z% = y\ kx 2 + (2y 2 + y) 2 . Under the non-linear change of variables 



yi 



W+ 371/1 + (37 -2) + 



2 ,4-3 7 



[x + y ) + -y 2 y, 



2/2 -> 2/2 + 4t/iy 2 + ^l/i 2/ , 

X — > X + 2?/ 1 X, 

y ->y- 22/12/2 + 2y x y, 



(A.28) 



and keeping only terms up to second order, the system transforms to 



2/1 



-1/2- -2/11/ + 0(3), 



2/2=2/1- 2^2/ + C (3) , 
x = -xy + O (3) , 

y = 6( 7 -l)(2/? + 2/ 2 2 )-^x 2 -|: y 2 + 0(3). 



(A.29) 
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